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Abstract. We present the first set of results of solving the Hartree-Fock- 
Bogoliubov equations, which describe the self-consistent mean field theory 
with pairing interaction. Calculations for even-even nuclei are carried out on 
a two-dimensional axially symmetric lattice, in coordinate space. An impor- 
tant aspect of our method is the proper representation of the quasi-particle 
continuum wavefunctions, which are considered for energies up to 60 MeV. 
This stage is essential for a proper description of nuclei near the drip lines, 
due to the strong coupling between weakly bound states and the particle con- 
tinuum for such nuclei. High accuracy is achieved by representing the operators 
and wavefunctions using the technique of basis-splines. Calculations for Sn 
isotopes are presented to demonstrate the reliability of the method. 



1. Introduction 

The study of structures and reactions of nuclei far from stability has been one of 
the most active fields of nuclear physics in the past decade 0,|2). The microscopic 
description of such nuclei will lead to a better understanding of the interplay among 
the strong, Coulomb, and the weak interactions as well as the enhanced correlations 
present in these many-body systems. 

The experimental developments as well as recent advances in computational 
physics have sparked renewed interest in nuclear structure theory. In contrast to the 
well-understood behavior near the valley of stability, there are many open questions 
as we move towards the proton and neutron driplines and towards the limits in mass 
number. The neutron dripline represents mostly "terra incognita" . In these exotic 
regions of the nuclear chart, some of the topics of interest are the effective N-N 
interaction at large isospin, large pairing correlations and their density dependence, 
neutron halos/skins, and proton radioactivity. Specifically, we are interested in 
calculating ground state observables such as the total binding energy, charge radii, 
proton and neutron densities, separation energies for neutrons and protons, pairing 
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gaps, and potential energy surfaces. It is generally acknowledged that an accurate 
treatment of the pairing interaction is essential for describing exotic nuclei |^,|||. 

In Section 2 the general outline of the HFB formalism is shown in second quan- 
tization and its representation in coordinate space. Section 3 shows the reduction of 
the HFB equations for cylindrical coordinates. The last two sections include results 
and discussion of our HFB calculations. 



2. Hartree-Fock-Bogoliubov Formalism 

2.1. Basic outline of HFB equations 

The many-body Hamiltonian in occupation number representation has the form 
H = ^2 < i\ t \j > c\ cj + j < y'l i{2) \ mn > d l d j 6 ™ £ ™ + ■■• C 1 ) 



4 



Similar to the BCS theory one performs a canonical transformation to quasiparticle 
operators ft ft 



c 



(3 \ _ ( W V* 
ft )-\v T U T )\ gt 



(2) 

The HFB ground state is defined as the quasiparticle vacuum 

ft |$o > = . (3) 

The HFB ground state energy together with the constraint on the particle number 
is given by 

E(U) =< $o\H- AiV|$ > , (4) 

where 1Z is the generalized density matrix, function of the normal density p^.j and 
the pairing tensor Kij. We derive the equations of motion from the variational 
principle 

S [E(K) - tr A{K 2 - K)} = , (5) 
which results in the standard HFB equations 

[H,K] = , (6) 

with the generalized single-particle Hamiltonian 

where h and h denote the particle and pairing Hamiltonians, respectively and the 
Lagrange multiplier A is the Fermi energy. 
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2.2. HFB equations in coordinate space 

For certain types of effective interactions (e.g. Skyrme mean field and pairing 
delta-interactions) the "particle" Hamiltonian h and the "pairing" Hamiltonian h 
are diagonal in isospin space and local in position space, 

h(raq,r'a'q') = 6 q , q , S(r - r')^, (r) (8) 

and 

h(raq, v'a'q') = 6 q<q , S(r - r')h%, (r) . (9) 

Inserting these into the above HFB equations results in a 4x4 structure in spin 
space: 



{h q - A) h q 



=eM~ ) , run 



with 



and 



Mt« WJ ' W T (r) h^r) 



(12) 



Because of the structural similarity between the Dirac equation and the HFB 
equation in coordinate space, we encounter here similar computational challenges: 
for example, the spectrum of quasiparticle energies E is unbounded from above 
and below. The spectrum is discrete for \E\ < —A and continuous for \E\ > —A. 
For even-even nuclei it is customary to solve the HFB equations with a positive 
quasiparticle energy spectrum +E a || and consider all negative energy states as 
occupied in the HFB ground state. 



3. 2-D Reduction for Axially Symmetric Systems 

For simplicity, we assume that the HFB quasi-particlc Hamiltonian is invariant 
under rotations R z around the z-axis, i.e. [H, R z ] = 0. Due to the axial symmetry 
of the problem, it is advantageous to introduce cylindrical coordinates (<p,r,z). 

It is possible to construct simultaneous eigenfunctions of the generalized Hamil- 
tonian TL and the z-component of the angular momentum, j z 

jz ipn,n, q {<i>,r,z) = htt ipn,si,q(<l>,r,z) , (13) 
with the quantum numbers il = ±^,±|,±|,.... The simultaneous quasiparticle 
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eigenfunctions take the form 



*Pn,n,q{<j),r, z) 



,r, z) 



/2tt 



/ e m-^ ^ nq{r ,z,T) \ 

\ ^ ] U%, q {r,z,i) I 



(14) 



We introduce the following useful notation 

U^ q \r,z) = ^] q {r,zA) , L^ q \r,z) = ^] q (r,z,i) 



(15) 



From the vanishing commutator, \H,j z ], we can determine the (^-dependence 
of the HFB quasi-particlc Hamiltonian and arrive at the following structure for the 
Hamiltonian 



h{(f>,r,z) = 
and the pairing Hamiltonian 

h(4>,r,z) = 



fr' TT (r,z) 
h' n (r, z) 



a-** h' n (r, z) 
h 'll ( r ' z ) 



h' n (r, z) e-** h' n (r, z) 
+** h' (r,z) h> (r,z) 



(16) 



(17) 



Inserting equations ( |l6| ) and ( |l7j ) into the eigenvalue Eq. ( |l0| ) , we arrive at the 
"reduced 2-D problem" Q in cylindrical coordinates: 



f(h' n -\) 


TJ. 


h'n 


h'n 


h 'n 


(^U-A) 


h'n 


h'n 






-(h' n -X) 


— h' 

n n 




h'll 


-h'n 


'(h'n ~ 



A)/ 



(U (1) \ 

r(l) 
n n,Q,q 

U n,Cl,q 

\ T {2) I 



E, 



(u (1) \ 

' U n,n,q I 
Al) 

n,U,,q 
\ r(2) / 



Here /i', h! , [/'s, and L's are all functions of (r, z) only. For a given angular mo- 
mentum projection quantum number Q, we solve the eigenvalue problem to obtain 
energy eigenvalues -En,fi, 9 and eigenvectors Vyn, <? for the HFB quasi-particle states. 

3.1. HFB Hamiltonian using the Skyrme interaction 

Using the Skyrme form for the particle Hamiltonian, h q , we can write 



V — V + U q + Uc- iB q - (V x a) 

2m* 



where the effective mass is defined by 

h 2 



h 2 



2m* 2m q ' 4 



+ - A (t l +t 2 ) P +\(t 2 -t 1 ) Pq 



(18) 



(19) 
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Applying the cylindrical form of the Laplacian operator to the standard form of the 
wavefunction in Eq.(^), and invoking the axial symmetry of / we find 



hi o 
o t 22 



where the elements are given by 



(20) 



t - /( — + -—- 

\ dr 2 r dr 



t 2 2 = f 



d 2 ld_ 

dr 2 r dr 



(fi - 1/2) 



(0 + 1/2) 



d^_ 

dz 2 

&2 
dz 2 



df_d_ + d£d_ 

dr dr dz dz 

df_d_ d/_d 

dr dr dz dz 



(21) 
(22) 



/ being the effective mass. The local potential terms could also be cast into a 
matrix form 

' wii 

V 2 2 



(23) 



where 



«n = v 2 2 = U q + U c 



(24) 



The individual terms are constructed simply by summing appropriately weighted 
scalars as indicated by Eqs. (^) and including the Coulomb potential and the Slater 
exchange (26) term: 



U q = t 



1 



1 + o^o ) P - o + x o ) Pq 



1 



-(t 1 +t 2 )T--(t 1 -t 2 )T q 



1 fl 



ap 



a-l \ " 2 i r, a 



1 



t 4 [V ■ J + V • Jg] 



2 

+ g(«2-3t 1 )V a p + ^(t 2 + 3t 1 )V a p, 



j3„/ Pp( r> ) 2 



J \T-V 

The Hartree-Fock spin-orbit operator 



1/3 



\Pp(*)\ 



1/3 



iB q - (V x a) 



(25) 



(26) 
(27) 
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could similarly be cast into the form 



with 



Wn W12 
W 2 1 W22 



= B, 



n - 1/2 



Wu = 




+ 1/2 
r 


Or 


^4 

oz 


W21 = 




-1/2 
r 


+ B°- 
or 


oz 


W22 = - 




1/2 

1 







(28) 



where 



B r — B q • e r = V r y (/? + p q ) 



B z = B cl -e z = X7 z t -±{p + Pq ) . (29) 



3.2. Densities and currents 

While the form of the single particle Hamiltonian remains the same as Skyrme HF 
Hamiltonian, the densities and currents need to be written in terms of the quasi- 
particle wavefunctions. We obtain the following expressions for the normal density 
p q (r) and pairing density p q (r), which are defined as the spin-averaged diagonal 
elements 

Pi( r ) = J2p( rcr< l,rcrq) = foA^Q) ^(rcrg) (30) 



and 



Pi 



( r ) = p( r ^ r<j <i) = -J2J2 ( t )2 ^ raqS > ^*A r<j i) 



(31) 



The physical interpretation of p q has been discussed in the quantity [p q (r) AV/2] 2 
gives the probability to find a correlated pair of nucleons with opposite spin projec- 
tion in the volume element AV. 

Using the structure of the bi-spinor wavefunctions defined earlier we find the 
following expressions for the normal and pairing densities. 



p q (r,z) = i- 2 £f ] fj [\U% q (r,z)\ 2 + \L^ q (r,z)f 



n>o/ E n >o 



(32) 



1 / f2 mflX \ E inax 

\ n>o/ E n >o 



U%(r, z)U$£(r, z) + L%{r, z)L%{r, z) 



(33) 
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Similarly, the kinetic energy density and the divergence of the current become 

T,(r) = V-V7> g (r,r')| r=r , = £ £ | V fa, a {raq)\ 2 , (34) 



2tt 



2 £, . 

\ n>o / e„>o 



(n- 1/2) 5 



u. 



(2) 



+ 



(ft + 1/2) 2 



dr 



dL 



(2) 



dr 



dU : 



(2) 



32 



dL. 



.(2) 
J n£lq 



(2) 



<9z 



(35) 



V • J,(r) 



2tt 



2 E E 2 

\ I2>0 / E n >0 



filing ^Ljg ^isjg filing 



dr dz 



dr 



dz 



+ 



» -1/2 (2) f < ^in/ 



dr 



dz 



n + l/2 T (2) ( dU ua q , 

'• . ii.. I ~ r 



(2) aL (2) 



J ni1q 



dz 



dr 



The total number of protons and neutrons is obtained by integrating over their 
densities 



err p 9 (r) — 2n rdr 
o 



p 9 (r, z) 



(36) 



Finally, we state the normalization condition for the four-spinor quasiparticle wave- 
functions as 

(37) 



J d 3 r V&i,( r ) ^( r ) 



which leads to 



|E#LM| 2 + l4lM| 2 + A(r,z)\ 2 + \L%(r,z)\ 



1 . 



3.3. Pairing interaction. 



If one assumes that the effective interaction v p J ir is local, 

.,(2) 



<al(^,r' - ^jnVi.ra' - 4) = 5(n' - r) ^ <5(r 2 ' - r') ^^(ra, r' - o 1 ) 



one finds the following expression for the pairing mean field 
ft(rc7,rV) = V p (rcr,r' - er') p(rcr,rV) . 

For the local pairing interaction V p we utilize 

y p (ra, r' - a') = V 6(r - r') 6 a>(T , F{r) . 



(38) 
(39) 
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This parameterization describes two primary pairing forces: a pure delta interaction 
(F = 1) that gives rise to volume pairing, and a density dependent delta interaction 
(DDDI) that gives rise to surface pairing. The DDDI interaction generates the 
following pairing mean field for the two isospin orientations q = ±s 



1 



h q (r*,r'*') = ^-V^p q {v)F(v)5(v-v')5^, . 



The pairing contribution to the nuclear binding energy is 



P _ Tp(P) 17' 

Impair — ^p a i r T ^ p air 



(fir 



(P) 



~ 2 
-Pn 



Vr 



(™) 



F(r) 



(40) 



(41) 



3.4. Lattice representation of spinor wavefunctions and Hamiltonian 

For a given angular momentum projection quantum number f2, we solve the eigen- 
value problem on a 2-D grid (r a , zp) where a = 1, N r and i3 = 1, N z . The four 
components of the spinor wavefunction ip(r, z) are represented on the 2-D lattice by 
a product of Basis Spline functions Bi(x) evaluated at the lattice support points. 
Further details are given in Ref . Oj . 

For the lattice representation of the Hamiltonian, we use a hybrid method 
p0|,p^| in which derivative operators are constructed using the Galerkin method; 
this amounts to a global error reduction. Local potentials are represented by the B- 
Splinc collocation method (local error reduction) . The lattice representation trans- 
forms the differential operator equation into a matrix form 



N 



(/!= 1,...,JV) , 



(42) 



with N = 4 x N r x N z . The method of direct diagonalization with LAPACK is 
implemented to solve this eigenvalue problem. Our HFB code is written in Fortran 
90 and makes extensive use of new data concepts, dynamic memory allocation 
and pointer variables. The code uses as a starting point the result of a HF+BCS 
calculation, which makes HFB converge substantially faster. 

Since the problem is self-consistent we use an iterative method for the solution. 
The Fermi level, A, is calculated in every iteration by means of a simple root search, 
and used for the next iteration. This process is done until a suitable convergence 
is achieved. The quasiparticle energies and corresponding wavefunctions are calcu- 
lated up to 60 MeV. In practice a cutoff at this energy is imposed, but this limit 
can be set higher if necessary. Further details will be given elsewhere 1 16 . 



4. Results 

In table [l] we diplay the results of calculations for two tin isotopes 120 Sn and 150 Sn. 
In the calculations of 150 Sn we used a box size (Rbox) of 20 fm. The numerical 
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mesh includes 17 and 34 points in r and z direction respectively. These points are 
geometrically distributed, giving more data points in the central region, where the 
particle and pairing densities are denser. The maximum SI number was for this 
case. For the 1-D calculations, Rbox = 30/m, linear spacing of 0.25 fm and j max 
of 4r was used. 

Table 1. Comparison of calculations for tin isotopes. The 1-D calculations were 
made by Dobaczewski ||J using the SkP interaction for 120 Sn, and SkM* for 150 Sn. 
Calcu 





120 Sn 


150 


Sn 


Observables 


Exp. 


1-D 


2-D 


1-D 


2-D 


Binding Energy (Mev) 


-1020 




-1021 


-1155 


-1157 


Fermi Energy (n) (Mev) 




-7.94 


-8.39 


-1.86 


-1.99 


Fermi Energy (p) (Mev) 






-7.58 


-16.13 


-15.95 


Pairing Gap (n) (Mev) 


1.378 


1.256 


1.85 


1.55 


2.10 


Pairing Gap (p) (Mev) 




0.00 


0.00 


0.00 


0.00 


Beta 2 


-0.004 




-0.008 




-0.001 



According to the results shown in table |l| for Sn isotopes, the agreement is 
evident with respect to the 1-D calculations by Dobaczewski, although 120 Sn cal- 
culations were done with different Skyrme forces. We encounter some differences 
in the pairing gap. This is likely to be explained by the dependence of the pair- 
ing gap on the box size, which was 30 fm for 1-D calculations but 20 fm for 2-D 
calculations. 



Total-Neutron-Potential 



Neutron-Pairing-Density 




Fig. 1. Neutron potential and pairing density for 50 Sn 



In figure [l] we show the total neutron potential and the neutron pairing density. 
As we can see the pairing effect in the case of 150 S'n is appreciable. 
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5. Conclusions 

We have seen that our calculations are close to the experimental values, and that 
our two-dimensional, axially symmetric code agrees with the calculations of the 
one-dimensional HFB code Q for spherical nuclei. 

The axial symmetry imposed in our HFB code is expected to be well suited in 
describing deformed nuclei far from stability. Our HFB approach for such nuclei 
works especially well in treating strong coupling to the continuum, which was shown 
to be crucial for obtaining convergence. 
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